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Abstract 

The parameter dependence of the various attractive solutions of the three variable nonlinear 
Lorenz model equations for thermal convection in Rayleigh-Benard flow is studied. Its bifurcation 
structure has commonly been investigated as a function of r, the normalized Rayleigh number, at 
fixed Prandtl number a. The present work extends the analysis to the entire (r, a) parameter plane. 
An onion like periodic pattern is found which is due to the alternating stability of symmetric and 
non-symmetric periodic orbits. This periodic pattern is explained by considering non-trivial limits 
of large r and a. In addition to the limit which was previously analyzed by Sparrow, we identify 
two more distinct asymptotic regimes in which either ajr or a 2 /r is constant. In both limits the 
dynamics is approximately described by Airy functions whence the periodicity in parameter space 
can be calculated analytically. Furthermore, some observations about sequences of bifurcations 
and coexistence of attractors, periodic as well as chaotic, are reported. 
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I. 



INTRODUCTION 



The nonlinear coupled set of three ODEs 



x = —ax + ay , 



y = rx — y — xz , 



z = —bz + xy , 



(1) 



was first derived by Berry Saltzman and Edward Lorenz |2j as a model for thermal 
convection between two plates perpendicular to the direction of the earth's gravitational 
force. It has three control parameters. First, there is the normalized temperature difference 
r = AT/AT C = Ra/Ra c between the hotter lower plate and the colder upper one. AT C is the 
temperature difference at the onset of convection. Ra = ga p L 3 AT/ uk, the Rayleigh number, 
is the dimensionless temperature difference, g denotes the gravitational acceleration, a p the 
isobaric thermal expansion coefficient, L the distance between the plates in upward direction, 
and v, k characterize the kinematic viscosity and thermal diffusivity, respectively. 

Second, there is the Prandtl number a = v/k, describing the ratio of the two molecular 
time scales r v = L 2 /u and r K = L 2 / n for momentum and energy transport in the system. 
a also measures the ratio of the widths X u and At of the kinematic and the thermal boundary 
layers. Large Prandtl numbers a indicate a broad kinematic boundary layer containing 
a small thermal one, fast molecular momentum exchange and small convective Reynolds 
numbers. Small Prandtl numbers characterize fast molecular thermal exchange, a broad 
thermal boundary layer containing a thinner kinematic one and larger convective Reynolds 
numbers. Typical values for a are 0.7 in air (at room temperature), N 2 , O2, He, and other 
gases, about 7 in water (room temperature), 600 through 3000 in typical organic liquids as 
e.g. dipropylene glycol (600), or a ~ 7250 in glycerine. Small Prandtl numbers are realized 
in mercury (a « 0.025) or liquid sodium with a m 0.005. 

The third parameter in eqs. is the coupling strength b, originating from the nonlin- 
earity (u ■ V)T of the Boussinesq equation. It is of order 1 and conventionally chosen as 
b = 8/3; we adopt this value here. 

The variables x(t),y(t), and z(t) describe, respectively, the amplitude of the velocity 
mode (the only one left in the model), the amplitude of the corresponding temperature 
mode, and the amplitude of the mode measuring the nondimensional heat flux Nu, the 
Nusselt number, x, y characterize the roll pattern, the heat flux is 



Nu = 1 + 2z/r . 



(2) 
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For < r < 1 there exists one stable fixed point Po (xq = yo = Zq = 0), with Nu = 1, 
corresponding to pure molecular heat conduction. For r exceeding 1, convective rolls develop, 
described by two stable fixed points 

P± = (x±, y±, z±) = (±y/b(r-l), ±y/b(r - 1), r - l) (3) 
hence Nu = 3 — 2/r. The trajectories in phase space T = {x, y, z} are spirals towards these 



fixed points. Stability is lost for 



r > r c = a . (4) 

a — b — 1 



For the canonical choice a = 10 it is r c = 24.7368.... At such large r the three mode approx- 
imation for the PDEs describing thermal convection has of course ceased to be physically 
realistic, but mathematically the model now starts to show its most fascinating properties, 
because a strange attractor develops for r > r c . With further increasing r, windows of regu- 
lar behavior emerge out of chaos where stable periodic orbits exist. These stability windows 
tend to become broader and broader until for sufficiently large r chaos seems to disappear 
and one type of stable periodic solution takes over. 

Of the numerous studies of the Lorenz model we cite as examples the monograph |3j and 
the more recent survey in Jackson's book |4|. To the best of our knowledge, the parameter 
dependence of the Lorenz model has mainly been analyzed with respect to its control pa- 
rameter r, the external forcing, and sometimes with varying b (as in the work which is 
cited in |4J). The Prandtl number a has usually been fixed to its canonical value a = 10, and 
as a rule, b has been taken as 8/3. The limit of large r at arbitrary fixed o was thoroughly 
studied by Robbins jf| and Sparrow ^j. 

In the present paper we study the parameter range where both r and a become large. 
The scheme of bifurcations is then a very different one. With a as large as 60 (Fig. one 
is just beginning to see a structure that unfolds on much larger scales (Fig. |2J). There is 
ongoing repetition of a few basic features, but these are different from the structure at low o. 
The line a s = Vbr divides these two regimes as we shall show. Our theoretical explanation 
of these findings focuses on a backbone structure of the dominant stable periodic attractors, 
but orbits with arbitrarily high periods, and chaos, are also observed as part of the pattern. 

An important feature of the Lorenz equations is their invariance with respect to the 
reflection 

S : x, y i— > —x, —y and z h- > z . (5) 



Evidently the z-axis is a S-invariant line. For r > 1 it is also the stable manifold of the 
fixed point Po- This characteristic symmetry property with respect to £ has implications 
on the nature of the period doubling scenarios which occur over and over again as the 
parameters r, a are varied: they do not follow the familiar scheme of the paradigmatic 
logistic iteration x h- > ax(l — x) but rather the somewhat more intricate scheme of the 
cubic iteration x t— > f(x) = (1 — c)x + cx 3 . This cubic mapping possesses the symmetry 
/ = S^/E, with £ : x i— > —x, and may be viewed as one building block for understanding 
the transition to chaos in the Lorenz system, see Hao et al. ?[ |s| for a detailed analysis of 
this aspect. 

The motivation of the present work was the simple question why the meteorologist Lorenz 
would have chosen the Prandtl number a = 10, which is adequate for water but less so for 
air with its order of magnitude a ~ 1. The reader will easily find that taking a this small, 
he will no longer find chaotic solutions. This is confirmed by a glance on eq. (J3J): with £7=1, 
the critical parameter r c = —5/2 lies in the unphysical range r < 0. If one chooses, instead, 
some cx > 6+1, implying r c > 0, the window structure in the r-dependence is similar to 
the case a = 10 although the details may be nontrivially different. But one hardly sees how 
these differences in the r-behavior for some discrete a-values are connected. This motivated 
us to vary a in sufficiently small steps and study the 2-dimensional (r, cr)-space of control 
parameters in a larger domain than has hitherto been done. However, quite apart from such 
a physical motivation, the mathematical complexity of this innocent looking set of equations 
warrants an investigation of their asymptotic behavior where it has not yet been explored. 

In Section |H] we present our parameter survey (phase diagrams) and discuss its main 
features. Then, in Section IIIII the nature of two kinds of periodic orbits is described whose 
stability regions define the dominant regular features of the large scale phase diagram; an 
analytic derivation of these observations is given in Section HVl More complex aspects of the 
Lorenz dynamics such as bifurcations of periodic orbits, chaos, and coexistence of attractors, 
are discussed in Section El We point out regularities, and similarities to other well known 
systems, but do not give detailed explanations of these features. 
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II. TWO-PARAMETER ANALYSIS 



A glance at Figs. ^ and Ogives an impression of the bifurcation scheme on large and very 
large scales of r (abscissa) and a. The parameter b, as mentioned already, is kept at its 
standard value 8/3. Colors indicate the existence of attractors of a certain type, which does 
not preclude the coexistence of others. More precisely, the figures are parameter scans where 
each point (r, a) is colored according to which type of attractor is approached from a fixed 
initial point (xo, yo, zo). Such pictures are called final state diagrams, or phase diagrams. 

The simplest attractors are the fixed points P±, see 0. Their stability is lost along the 
curve r = r c (a), see (@J), or rather by its inverse function a = <J c [r) with 

°c(r) = r -- h -^± l -^r*- 2(36 + 5)r + (b + 3) 2 . (6) 

For large r, the upper and lower asymptotes a u (r) and <Je{r) are, respectively, 

a u = r -2(6 + 2) and <r e = b+l. (7) 

The next order correction is =f(6+ l)(26 + 4)/r, describing the approach of the upper asymp- 
tote from below and the lower one from above. 

The (r, <r)-region where the points P± are stable according to this linear analysis is colored 
black, the rest is left white, before the other colors are added in Figs. and |21 These colors 
code for the periods of attractors which are determined by the procedure described below: 
yellow for orbits of period 1, red for period 2, green for period 4, blue for period 8, grey for 
period 3, olive for period 6. 

For each (r, a) on a grid of 610 x 620 points, an orbit was started with initial conditions 
(x ,y ,z ) = (0.001, 0.001, r — 1). A transient regime of 80 000 iterations was allowed to 
elapse before it was assumed that an attractor was reached. The trajectory was then ana- 
lyzed in terms of its intersections with the (x, ?/)-plane at z = r — 1 which is an appropriate 
Poincare surface of section. The period p of a periodic orbit is defined as half the number 
of intersection points of the orbit with that plane (in one half the trajectory intersects from 
above, in the other half from below). Fig. El gives two examples for illustration. To deter- 
mine the period we used a rather simple scheme. A disk of radius 10 -3 was placed around 
each point of intersection unless it fell in a disk already generated. When the number of 
patches so obtained remained finite, like in the left part of Fig. El the orbit was assumed to 
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FIG. 1: (?", d)-phase diagram in the range < r < 409 and < a < 60. Colors indicate the type 
of attractor reached from the initial condition (xq, yo, Zq] = (0.001, 0.001, r — 1): fixed points P± 
(black), orbits of period 2 (red), period 4 (green), period 8 (blue), chaotic attractors (white). This 
picture is a zoom into the lower left corner of Fig. [2J The line a s := Vbr in color magenta marks 
a transition range between different global structures for smaller or larger a. Further details are 
discussed in the text, Section IV Al 

be periodic. When, on the other hand, the number of disks seemed to grow indefinitely, like 
in the right part of the figure, the orbit was considered to be chaotic. In practice we waited 
for more than 100 intersections, but only periods up to 8 form clearly identifiable bands on 
the scales of the figures. When a periodic orbit was detected, the color of the grid point 
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FIG. 2: The parameter ranges are < r, a < 2000. The color code is the same as in Fig. ^ 
Yellow bands indicate stability of orbits with period 1 (not present in the zoomed lower left corner, 
cf. Fig. [J). The dots in the red and yellow bands along the line a = 0.75 r are, from bottom left to 
top right, Si, Ai, S2, A2, . . . , defined according to Q44|) for rjq{m = 0.75). With their upper parts, 
the bands penetrate well into the range above a u = r — 28/3, cf. ((7J) (the line can be identified as 
the border to the black region through the bands), indicating coexistence of periodic orbits with 
the stable fixed points P±. 
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FIG. 3: (Color online), (al): (y, z)-projection of an asymmetric attractor of period 3 (abscissa y, 



ordinate z). r = 1744, a = 110, initial conditions xq = yo = 10 
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1. The fat points are 



the fixed points P±. (a2): The same for a chaotic attractor; r = 1744.9, a = 115.2. (bl) and (b2): 
Corresponding Poincare sections with the (x,y)-plane at z = r — 1; intersection points are grey 
(green) if from below, black (red) if from above (abscissa x, ordinate y). 



(r, a) was overwritten with the color of the code given above. Otherwise the original black 
or white color was left unchanged. Hence chaotic attractors appear as black or white points. 

The method produces two artifacts. First, it cannot identify the fixed point attractors P±; 
hence, in the black region it does not distinguish between fixed points and chaotic attractors. 
Second, when the relaxation towards an attractor is very slow, the 80 000 transient iterations 
may not be sufficient to reach it, and the method produces more than 100 patches. Hence 
this situation is not distinguished from chaos. It is observed in connection with critical 
slowing down in the neighborhood of a bifurcation of the attractor. Examples are the white 
bands in the red regions at bottom and top of Fig. the narrow seams between red and 
green regions, and the lines inside the two green bands in the chaotic range. The nature of 
the corresponding bifurcations will be elucidated in the following section. 

Since only one initial condition is used, the procedure does not determine whether dif- 
ferent initial conditions would lead to different attractors for the chosen parameters (r, a). 
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The question of coexistence of attractors will be addressed later, in Sec. IV Bl At this point 
we can assert already that where a colored band extends into the black region, like in the 
lower right and upper left parts of Fig. ^ the corresponding periodic orbit coexists with the 
fixed points P±. 

The arrangement of bands in Fig. Q exhibits a qualitative structural symmetry with 
respect to a line which has been added in magenta color. This line, mathematically defined 
as cr s = Vbr, separates a large-cr structure from a kind of mirror image in the low-cr part of 
the phase diagram. It is conspicuous how every band below cr s (r) seems to have a counterpart 
on the other, upper side. But, as a glance on Fig. El shows, this "symmetry" disappears on 
a larger scale: there is no repetitive structure in the low-cx region, in contrast to the large-a 
range. 

No orbits of period 1 (yellow bands) are observed in the (r, <r)-range covered in Fig. d 
However, in Fig. El which covers a much larger (r, <j)-range, an onion like repetitive pattern 
of yellow (period 1) and red (period 2) bands dominates in the upper part. These bands 
are separated from each other by chaotic regions which also contain small bands of higher 
periods. The bands (above a s ) have an hyperbolic shape whose upper branches seem to have 
asymptotic slopes of approximately 1.3 while the lower branches will be shown in Section ITVl 
to have a limiting behavior a oc Vbr like the "bottom line" <r s (r). The upper branches 
penetrate significantly into the region a > o u where the fixed points P± are stable, see 
eq. ((Zj). This means there is coexistence of fixed point attractors and periodic attractors. To 
generate Fig. El the initial conditions were xq = yo = 10" 5 and z = r + 1000. Choosing the 
initial conditions as in Fig. the corresponding picture (not shown here) looks different in 
that the greater part of the white bands of chaos is replaced by colors yellow or red. This 
reflects the coexistence of chaotic and periodic attractors. 

The bands have infinite vertical slopes along the line a v (r) ~ 0.75r. On this line, we 
have marked points Si, Ai, S2, A 2 , S3, ... in the red (period 2) and yellow (period 1) bands, 
respectively. The names of these points will also be taken to denote the corresponding bands 
because the attractors of a given band are topologically identical. Their nature will first 
be described phenomenologically in the subsequent section, where we shall also see how it 
changes from one band to the next. An analytic explanation of the findings will then be 
given in Section ITVl 
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III. ATTRACTIVE PERIODIC ORBITS 



In this section we explore the nature of the dominant periodic attractors, i. e., of the orbits 
with periods 1 or 2, by means of numerical and graphical investigation. Their bifurcations 
into orbits of higher periods, as well as into strange attractors, will be discussed in Section IV! 

To begin with, consider the attractors corresponding to the (r, <r)-values of A3 in Fig. 0] 
and of S4 in Fig. in the three standard projections. Comparison with the Si, A 1; S2, and 
A 2 attractors (see Figs. ITIandlTBj) shows that on the scales used, — 0.4r < x,y, z — r < 0.4r, 
all A n attractors look alike, and so do the orbits of S n . Furthermore, the attractor of the 
parameter point S n (for given n) looks like a combination of the asymmetric orbits of A n 
and its S-mirror image SA n . 

The time course of the variables x, y, and z is plotted in Fig. El both for the A3 attractor 
(top) and the S 4 one (bottom). 

The combined picture shows the typical behavior of a relaxation oscillator. Two different 
regimes can clearly be distinguished: First, the variable z decreases slowly from an upper 
value r(l + p) to a lower value r(l — p), with p « 0.35, while x and y are virtually zero. 
Then, in a very fast process, x and y start out on a large excursion with amplitude rp, while 
z jumps back to its upper value. The (y, z)-projection exhibits an almost perfect half circle; 
the (x, ^-projection shows that during the excursion, x lags somewhat behind y. 

From Figs. HHHl it is not obvious what happens during the slow phase where the orbits 
move down along the z-axis. This question is clarified by the sequence of (x, ^-projections 
in Figs. [7| and |SJ The individual pictures in Fig. [7|show the attractors of Si, Ai, S 2 ; the other 
projections would be similar to those in Figs. 0] and EJ Taken together, the pictures make it 
clear that the orbits wind around the z- axis on spirals with decreasing R 2 = x 2 + y 2 . Coming 
from large values of R, near the maximum value of z, they enter a close neighborhood of 
the z-axis when they pass a point where y = 0. After a quarter turn in clockwise direction 
(looking down the z-axis) there comes a point where x — 0, at smaller values of z and R. 
What happens next depends on the index n of S n or A n . In the case of Si the orbit quickly 
leaves the (x, |/)-origin, cf. the inset in the left part of Fig. [7[ it performs a large loop and 
returns to the (x, y)-origin region from the other side, but again near the maximum of z. 
All in all it encircles the z-axis once, its winding number W with respect to the z-axis 
being equal to the number N — 1 of zeros in y (and correspondingly in x) upon each close 
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FIG. 4: Stable asymmetric periodic orbit corresponding to the parameters (r, a) = (1322, 992) 
of the point A3 in Fig. [21 From left to right, the attractor is shown in (x,y), (x,z), and (y, z) 
projection (the abscissa is always mentioned as first coordinate, the ordinate as second). The dots 
mark the projections of the fixed points P±, see (jSJ). The arrows indicate the direction of flow; 
the decrease of z takes place very close to the z-axis. The attractor coexists with its S-partner 
which is obtained by reflection of x and y. Together these two attractors look like the symmetric 
attractor in Fig. El The scales are —530 < x,y, z — r < 530. 




FIG. 5: Stable symmetric periodic orbit corresponding to (r,a) = (1552,1164), the point S4 in 
Fig. [21 Same projections as in Fig. [IJ The scales are —620 < x,y, z — r < 620. 

encounter. In the case of Ai, the orbit stays longer near the z-axis, with a second zero of y 
and also of x, each with smaller values of R (and z); this is illustrated by the sequence of 
zooms in Fig. |HJ Its winding number with respect to the z-axis is also W = 1, but now there 
are two of these orbits, shown in different colors. To keep track of the winding number W 
for increasing band label n, we may focus on the number N of zeros of y (or x) upon close 
encounter with the z-axis. Detailed observation shows that this number is 



W(S n ) = N(S n ) = 2n - 1 and W(A n ) = ±JV(A„) = n. 



(8) 



Fig. IH1 shows a logarithmic plot of |x| and \y\ for A 3 (N = 6, W = 3) and S 4 (N = 7, W = 7). 
We observe three different regimes of time dependence. Soon after z reaches its peak value, 
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FIG. 6: (Color online). Time development of x, y, and z within one unit of time t (abscissa), 
for the attractors of A3 (top) and of S4 (bottom). The upper parts show the z-coordinates, the 
lower parts x and y which are nearly indistinguishable on this time scale. The vertical scale is 
— 0.4r < x,y, z — r < 0.4 r. 




FIG. 7: (Color online), (x, y)-projections of (from left to right) the attractors of Si, Ai, S2. In the 
case of Ai the two asymmetric orbits are shown in different color. The insert in the left picture is 
a blow-up of the neighborhood of the origin. The (r,a) values are Si = (172,130), Ai = (402,302), 
S 2 = (632,474); the scales are -0.4 r < x,y < 0.4 r. 
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FIG. 8: (Color online). Successive blow-ups of the box in the middle picture of Fig. [7| The 
sequence demonstrates that the two orbits of Ai do not reach the origin; rather each of them winds 
once around the z-axis. 

at a time to, the first zeros of y and x occur, visible as logarithmic singularities. The first 
regime of oscillatory exponential decay of x and y ends at a time t\. In the second phase 
there are no further x- or y-zeros, but the z-decay continues with deceleration up to a time 
t 2 where the distance R of the orbit from the z-axis reaches a minimum. The third regime 
is one of explosive, i. e., faster than exponential increase to the original values of x, y and z. 

The spiraling behavior in the first regime seems reminiscent of "Shilnikov Chaos" j^]]. 
However, in this case all the relevant dynamics stays far away from the fixed point. The 
invariant line x = y = is part of the stable manifold of the origin, but this is not important 
for our arguments. The main ingredient is that it is attractive for z > r (i. e. up to time 
t-i) and repulsive for < z < r. Long before the orbit gets close to the origin it leaves the 
z-axis along its unstable eigenvector, with x — y ~ 0. 

The global structure of the (r, a) parameter space, as exhibited in Fig. is obviously 
related to the existence and stability of such orbits, with winding numbers W(A n ), W(S n ) 
increasing as r and a (and thereby n) grow. The observations to be explained in the 
subsequent theoretical section are the following: 

1. Repetition of bands: Along the line a = mr with a fixed ratio m, there exists a period q 
in r (and consequently mq in a) such that if an attractor of type A n (or S n ) exists 
for values (r n ,cx n ), then a corresponding attractor of type A n+N (or S n+ Ar) exists for 
(r n + Nq,a n + Nmq). 

2. m-dependence of r-periods q: From the shape of the bands it may be inferred that 
q(m) has a minimum near m = 0.75. This means, the periodic repetition of the bands 
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FIG. 9: (Color online). Time development of log|x| (red), log|y| (blue), and above them z in 
grey (green) within 0.4 time units, for the attractors of A3 (top) and of S4 (bottom). The vertical 
scales are the same as in Fig. El x and y are plotted as 30 log \x, y\ + 600. Thus the x- and y-zeros 
appear as logarithmic singularities. 

is fastest (in r) along the line a = 0.75r. 

3. Size scaling of dominant orbits with r: The size of these orbits, i.e., the radius p 
of their (y, ^-projections measured in units of r, is independent of r: p « const for 
fixed m. 



Thus far we have not addressed the question of how p depends on m, so let us see next what 
theory has to tell; the result is presented in Fig. ^2 where p = p/m is plotted as a function 
of m. Of course, the above statements only hold for large r and a, in practice for r > 100 
and a in the bands above a, = Vbr. 



14 



IV. ANALYSIS OF THE PARAMETER SPACE STRUCTURE 



We begin with a discussion of appropriate scalings of the entire set of Lorenz equations JTJ) 
when both parameters r and a become large. Their asymptotic behavior with r — > oo de- 
pends on how a is assumed to change simultaneously. Robbins in (| and also Sparrow 
considered the case a = const. Here we study the behavior of the Lorenz equations with 
Prandtl numbers increasing as a oc r s , with some fixed exponent s > to guarantee increas- 
ing a, but s < 1 for o to stay below the upper stability border a u (r) = r — 28/3. 

If r and a become very large, the motion in x and y is expected to become very fast 
because d/dt in the equations of motion becomes large. Thus we introduce a rescaled, 
stretched time r by defining dr = e _1 dt with a scaling factor e~ l oc r e , e > 0. 

Let us now also rescale the variables x and y which measure, respectively, the amplitudes 
of the role velocity and temperature. As for the rescaling of z a glance on the second Lorenz 
equation shows that y is driven by (r — z). This fact is reflected in the amplitudes z of 
the numerical trajectories, or in the behavior of the fixed points P±, both as functions of r. 
Thus we define 

x = a£, y = Pv, z-r = 7C (9) 

The Lorenz equations then read, with £' = d£/dr etc., 

3 ex. T ct3 

£' = -ea£ + e*Z-ri, v ' = -eq - ey-C£, (' = -e6(C + -) + ^- (10) 
a p 7 7 

The freedom of choosing scaling parameters may now be used to normalize to unity the 

coefficients of the driving term in the first equation, and of the two coupling terms oc (£, and 

oc £77, which leads to 

eaB = a, e'ja = 3, sad = 7. (11) 

These equations allow us to express a, 3, and 7 in terms of a and e. We find first e 2 a 2 = 1, 
thus a = ±£ _1 , then 3 = ie^cr" 1 , and finally 7 = e~ 2 o~ x . The ± signs may be interpreted 
as corresponding to the symmetry E which is clearly preserved under the rescaling. We can 
therefore restrict ourselves to the plus sign and have 

a = e- 1 oc r e , (3 = 7 = e^a' 1 oc r 2e ~ s . (12) 

The asymptotic behavior will depend on how the exponent e compares to the other expo- 
nent 2e — s. There appear to be three distinct possibilities. One is that 3 and 7 grow faster 
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than a, i.e., 2e - s > e or e > s. The second is that all three variables x,y,z scale with 
the same factor, 2e — s = e or e = s. The third possibility e < s may be excluded with the 
argument that then the relaxation term — ea^ in the ^-equation would diverge. 

A representative case of the first possibility is s = and e = i,meaning a = const and 
e oc 1/y/r — > 0. This was the choice of Robbins and Sparrow |3]. All dissipative terms 
vanish in this limit, and the equations reduce to 

£' = 77, »/ = -«, C' = £7, (a = const, e = l/y/7). (13) 

This system has the two integrals A = |£ 2 — C an d B 2 = rf + £ 2 ; it is equivalent to the 
equations of a pendulum with length B and energy A. The values of the integrals A and B 
have to be determined from the e-correction to (J13j) . i. e., by considering their slow dynamics. 
In (21 0] it was shown how this can be done in terms of elliptic integrals. In that sense, this 
case has been solved. 

Let us therefore concentrate on the second possibility e — s. Then ea is independent of r, 
and we can take ea = 1. All scaling factors can now be expressed in terms of the Prandtl 
number a. The term — erj in the second Lorenz equation may be neglected as vanishingly 
small in comparison to with ££. Hence the rescaled eqs. fjTT)J) become 

b br 

= + V' = -CZ, C' = ^--C--, (e = sanda = r s ) . (14) 

a a z 

But what about the two last terms in the ^-equation? Their order of magnitude is br~ s and 
br l ~ 2s , respectively. Again, there are three possibilities: 

Case A. With the choice 1 — 2s = 0, or s = 1/2, the second term will tend to a constant and 
dominate the decreasing first term. We can write a = wy/br with some constant w. The set 
of equations becomes 

£' = — £ + 77, 7/ = — ££, C = £77 — m; -2 (e = s = \ and cr = wVbr) . (15) 

When s varies between 1/2 and 1, the same equations apply because the last term in the 
^-equation of ([Hjl dominates the second. Note that no comparison can be made between 
these two terms and £77. 

Case B. When s — 1, the two last terms in the ^-equation are of the same order of magni- 
tude r -1 . We write a = mr and obtain the set of equations 

£'=-£ + t?, 77' = -C£, C' = & ~ — (C + — ) (e = s = 1 and cr = mr) . (16) 

mr V m/ 
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Case C. When s > 1, a grows faster than r, so we leave the range where a < u u {r). Formally 
it is found that both terms besides £77 in the ^-equation may be neglected, but in contrast 
to eqs. (|13|) there remains a dissipative term in the ^-equation. The constant of motion 
B 2 = rj 2 + C, 2 constrains the solution to a cylinder. On this cylinder the equation reduces to 
the damped pendulum, and therefore almost all initial conditions decay towards the stable 
equilibrium on the respective cylinder. 

We thus have to discuss the asymptotic regimes A and B, and do this separately. To our 
knowledge, the two sets of eqs. (|15|) and (|T^|) have not been analyzed so far. As we shall 
see, their solutions lead to an understanding of the phase space structure displayed in Fig. El 
and of the numerical findings in Section 11111 



A. The limit of large r with a = w\fbr 

We start with the simpler case of the asymptotic Lorenz equations, the set of the rescaled 
eqs. Of all the complex dynamics that it contains, we only want to describe periodic 

orbits of the simplest kind. From our numerical studies we anticipate that these orbits have 
two clearly distinct phases. There is slow motion down along the z-axis, and a fast return 
to the original large value of z via a large loop away from the z-axis. 

The fast motion is readily described when we neglect the constant term in the ^-equation 
of Then, combining the equations for 77 and (, we obtain r]dr] = — (d( which after 

integration gives the circles 

V 2 + C 2 = P 2 (IT) 

with as yet undetermined radius p. A convenient parametrization of these circles is found 
with a rescaling of time, defining s as an arc length via ds = £(r)dr. The equations for rj(s) 
and ((s) are thereby linearized, with solutions 

77(3) = p sin s, C(s) = — p cos s . (18) 

This explains the observation that ( oscillates between two equal values of different sign, 
namely ±p. 

The equation for £, 

d^ = -^- r]) 01 d7 = -( 1 -^)J (19) 
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does not appear to be analytically solvable, but as long as £ is not too small, it describes a 
relaxation of £(s) towards rj(s). Hence, with a delay in r of the order of 1, £ follows 77 on its 
excursion to values of the order of p. 

Consider now the slow motion. It is given by the exact dynamics on the stable manifold 
£ = 7] = of the fixed point (0, 0, 0). The solution with initial point £(0) = p is £(t) = p — 
t/w 2 and reaches the lower value — p at time T = 2pw 2 . Linearizing the set of equations f)15|) 
around this motion, we obtain 

£' = -£ + 77, r/ = -(p-T/w 2 )£. (20) 

The r-dependent eigenvalues are — | ± a/A(t) with discriminant A = | — p + t/w 2 . Three 
r-regimes may be identified. Assuming p > 1/4, the discriminant is initially negative; the 
motion is oscillatory and decaying according to e~ T ^ 2 . This regime ends at the time T\ 
defined by A(ti) = 0, or t\ = w 2 (p — |). In the second regime, the motion continues to be 
decaying but is no longer oscillatory; this is the case between t\ and T2 = w 2 p = T/2 where 
C( r 2) = 0. The third regime, between T2 and T, sees £ and rj increasing at a growing rate. 

To proceed with the (£, ^-equations, we introduce a new variable X by separating the 
exponentially decaying part from £. With £(r) = e~ T ^ 2 X(r), the two linear equations of 
first order are transformed into 

X" - A(t)X = , or X" - -^1 = (21) 

w 2 

where f = r — T\, and X(f) = X{r — T\) lives on the interval —T\ < f < T — T\. This 
equation is the well known Airy equation. Its general solution is a superposition of the 
two Airy functions Ai(z) and Bi(z) where z = w~ 2 ^ 3 f. Both solutions oscillate for z < 0. 
For z > the solution Ai decays whereas Bi grows exponentially. We are interested in the 
solution X(f) = Bi(z), or 

£(f) = e- Tl/2 e^/ 2 Bi(w~ 2/3 T ) . (22) 

with —T\ < f < T — t±. It may be checked that 77(f) = e~ T l 2 Y{j) is governed by the same 
Airy equation and thus, up to a phase shift, has the same solution. An example of the 
function Bi is shown in the left part of Fig. where the parameters are those of S4 (the 
plots are made on the basis of eqs. ()40 p -(|42j ) ). The corresponding £(f) is shown in the right 
part of the picture, in logarithmic representation. 

18 



From here it is straightforward to determine the value of p, and the periodicity of the 
pattern of bands. 

We know from the analysis of the outer part of the attractor, see (j!8Jl . that the solution 
starts, for r = or f = — r 1; at about the same values £ and t] at which it ends for r = T: 

£(0) = £(T) . (23) 



Using the asymptotic behavior of Bi for large negative and positive arguments, 

(24) 



Bi(-N)« ^^cos(§|,r + J) 



V^\z\ 1/4 

1 

\ ■ • ■ 



Bi (M) « ^^7jexp{fH 3 / 2 } 



v^kl 1/4 

and the fact that the exponential terms dominate the algebraic factors by orders of magni- 
tude, we simply equate the exponents on the two sides of For f = —T\ eq. (|22J) gives 
£(0) = 0(1) (zero exponent); for f = T — T\, we add up the exponents in eqs. ()22|) and (|24|1. 
with z = w~ 2 / 3 (T — Ti), and obtain 

= - t - ^ + 1-" 273 ' 7 - ->i 3/2 - 4 r + 1> - -> 3/2 • < 25 » 

Inserting T = 2pu> 2 and T\ = w 2 (p — j) we obtain 

144p 2 = (1 + 4p) 3 . (26) 

The relevant solution is p = 1.35... =: p - 

To determine the values w = wn where the orbits have a certain number N of £- and 
//-zeros at each close encounter with the z-axis, we use the asymptotic behavior of the Airy 
function Bi(— \z\) as given in (|24|) . Consider first the symmetric orbits S n as they develop 
from time f = —t\ to T — T\. The coordinates £ and rj come in from large values (positive 
or negative) and go out to large values of the opposite sign. This means the phase shift 
between £(f = — 7i) and £(f = T — T\) is it; the cosine performs an odd number of half 
oscillations. For the asymmetric orbits A n , on the other hand, the phase shift is a multiple 
of 2n. Taken together, this leads to the "quantum condition" 

l^i 3/2 + f = iVvr, (27) 

where N = 1 corresponds to the first symmetric and N = 2 to the first asymmetric band. 
With p = po we get 



w N = 1.02V4JV - 1 , N = 1,2, ... . (28) 
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To wrap it up, the results for this regime of large r are the following. The iV-th band 
is centered around the line <jn = w^ybr] the size of the corresponding attracting periodic 
orbits in the original variables (x, y, z) is Po<Jn- But note this is only asymptotically true 
and can hardly be seen in the diagram of Fig. El except in the lower right corner where 
the Si-band is just bending over to follow the y/r behavior, r = 2000 gives the prediction 
o\ = 74.5v^3 = 129 which lies at the lower border of the Si-band. The bands corresponding 
to higher values of N have not yet entered this asymptotic regime. 



B. The limit of large r with a = mr 

Let us now study the set of equations (jTfij) 

£' = -£+77, v ' = -a, c = ^v-—U + - 

mr \ m 

where m = a jr is the constant parameter. Much of the analysis is similar as in the previous 
subsection, but the explicit occurrence of r in the equations adds complications. The general 
strategy is the same. For the outer, large £ and rj parts of the solution we neglect the term 
oc r -1 in the ^-equation and get identical results as in (fT7 |) -(fT3 |) . 

The slow motion along the z-axis is more complicated even though the ^-equation, ne- 
glecting the ^ry-term, decouples again from the rest: its solution is no longer linear in r. 
With the initial condition ({0) = p the approximate ^'-equation is solved by 

C(r) = (p+i)e-^-i, (29) 

and the time T needed to evolve from p to — p is 

m mr , 1 + mp 

T = — log- 1. (30) 

1 — mp 

Using this slow exponential decrease of C, (with the rate oc r _1 ) in the equation 77' = — C(t)£j 
we have again a time dependent linear system of equations for £ and m Its eigenvalues at 
fixed time are 

A=-i±VA with ^ = \-ar) = \ + i-{ P +i)e- {b,mr)T - (31) 

The discriminant A is negative in the r-range < r < T\ where T\ is defined by the condition 
A(n) = 0: 

mr 1 + mp 
b 1 + m/4 
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In this r-range the amplitudes £,77 oscillate and decay with damping rate — \. The second 
regime (further exponential decrease of £ and 77, but now without oscillations) ends at r 2 , 
defined by the condition that one of the eigenvalues be zero, i. e., C( r 2) = 0: 

7777* 

r 2 = — log(l + mp) . (33) 

From r 2 to T one eigenvalue is positive and the variables £ and 77 increase at a growing rate. 
This behavior is in complete qualitative agreement with Fig. El 

With the ansatz £(r) = e _T//2 X(r), the two linear equations of first order are trans- 
formed into X" — A(t)X = which is, in contrast to the previous case A, not the 
Airy equation because it is not linear in r. However, with the substitution r 1— > u via 
u = 2ry exp{ — (6/2mr)r} it may be reduced to the standard form of the Bessel equation. 
The £, 77 -solution is then a linear combination of the two Bessel functions of first and sec- 
ond kind of order — 2ry-^/l/4 + 1/m and argument u. For fixed, finite r but large r we 
have u ~ 2r^ — r, linear in r, but with a large additive constant. Thus in the limit of 
large r the asymptotics for Bessel functions with large negative order and large arguments 
is needed. The simplest way to obtain such asymptotics is to use an approximation by Airy 
functions from the very beginning when solving the A-equation. This is done by linearizing 
the discriminant A(r) in the original equation at an appropriate time r*, 

A(r) = A(r*) + £ (C(t*) + i) (r - r*) + . . . , (34) 

upon which we obtain the Airy equation for X(f) = X(r), 

X"-afX = 0, « = -^(C(r*) + -L), (35) 

where 

f = T -T* with t* = t* - A(r*)/a. (36) 

But what should we take for r*? A natural choice seems to be r* = n where A = and 
C = 1, whence 

a = «i = — ^-f 1 + ^) and tT — Ti (t*=ti). (37) 
tot V 4 / 

However, if we care to treat the decaying and the exploding part of the £, 77 -dynamics with 
equal weight, then it appears natural to take the minimum of £(r) = e~ T l 2 X{f) as the 
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reference point, i.e., r* = r 2 . Here C( r 2) = and A(t2) = \, so that from (J35j) and we 
obtain 

6 



a = a 2 



m 2 r 



and r* = r 2 



m 2 r 



r 1 + 0((m/4) 2 ) 



r' = r 2 j. 



(38) 



The last step involves r 2 = ri + r 2 — ri and an expansion of r 2 — t\ = ^ log(l + ^) to second 
order in the small quantity ^. The difference rj" — ri is itself of second order. Finally, when 
we are mainly interested in the initial phase of the decay, r* = ought to be a reasonable 
choice. Then C( r *) = P an d A(r*) = —(p— |) which leads to 

b 



a = «o 



m 2 r 



(1 + mp) and r* = n + 0((m/4) 2 ) (r* = 0) . 



(39) 



Again r* — ri is of second order in small quantities. Assuming that the m-corrections are 
not too large, we take r* = ri in all three cases and find the solution 



£(f) = e -^/ 2 e- f / 2 Bi(a 1 / 3 f) , 



(40) 



with — n < f < T — T\. An example of this time development is given in Fig. ITU1 where 
the parameters corresponding to S4 are taken. The left picture shows Bi(a 1//3 f) in the time 
range — n < f < 0, the right picture is a logarithmic plot of £(f) for an entire period. The 
qualitative agreement with the numerical integration of the Lorenz equations, see Fig. El is 
quite impressive. 




0.2 0.4 0.6 




FIG. 10: Plots of Bi(a Q /3 rt) (left) and log|f(t)| = log(exp(-mrt/2)|Bi(a2 /3 ^)l) (right), where 
t = f/r varies on the (r-independent) interval [— r±/r, (T—T\)/r\. The parameters are r = 1552, a = 
1164, corresponding to S4. 



As for the value of a, we take a 2 to determine p (a glance at Fig. El shows that the gross 
feature of the time development is nearly symmetric with respect to t 2 ), and ot\ or «o for the 
periodicity of the oscillations in the initial phase of the decay (which takes place between 
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t = and t = 7~i). This is why the a's are different in the left and right part of Fig. ITU1 
The right part shows that the periodicity condition £(0) = £(T), see ()23[) . is approximately 
fulfilled for S4 = (r,a) = (1552, 1164) if we take a = a%. The number of zeros, however, is 
not correct (6 instead of 7); had we taken a = cto, the values of m|£(t)| would have been 
different between beginning and end. On the other hand, the left part shows that with 
a = a>o we get the correct number of zeros (the same is true for a = a\). This reflects the 
fact that A(r) is not a linear function of r. 

Using the asymptotics of Bi, the periodicity condition £(0) = £(T) implies 

-ir+fv^(r-7i) 8 / a = o. (41) 

This leads to the following equation for the dependence of p on m, 

1+mp 4/ l + m/4\3/2 

V^log^ = - log- '—) . 42 

1 — mp 3 V 1 — mp / 

The left part of Fig. ^2 is a plot of the corresponding function p(m); the dots are derived 
from numerical test values of £ that were generated by integrating the full Lorenz equations. 
The agreement is satisfactory. In the limit of small m, obtained by expanding to lowest 
order in mp and m/4, we recover the result (|26[). With increasing m, the amplitude p of the 
(■-motion decreases from the value po = 1-35 to 0.493 at m = 0.75. For m of the order of 1 
the analysis becomes less and less reliable as the linearization of A(r), and hence the use of 
Airy functions, ceases to be a good approximation. 

When the value p(0.75) = 0.493 is compared to the numerically determined value of p, 
see e. g. Figs. 0] and we must remember that p was scaled with r whereas p is scaled 
with cr = mr; hence p = mp, and p(0.75) = 0.37 which agrees well with the numerical 
observations. 

To determine the periodicity of the (r, <r)-phase space pattern as we go along a line a = mr 
with fixed m in Fig. El the same arguments as in the derivation of (|2T|) lead to 

§V^i /2 + f = Ntt. (43) 

But again the question arises as to which a should be taken. Since p was calculated with a 2 
it seems consistent to do the same here; the result is the lower of the three curves in the 
right part of Fig. ^2 It might be argued, however, that the behavior of the solution between 
t = and Ti is rather determined by the slope of A(r) in that range, so we took a.\ and a 
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FIG. 11: Left: The amplitude p(m) = p{m)/m of the ^-motion as given by (|42|). The points 
are obtained by numerical integration of the full Lorenz equations. Right: The frequency uj(m) = 
2ir/P(m) of the repetitive pattern in the parameter r, as a function of m. From top to bottom, 
the curves axe calculated with. c n = c , ci, and c 2 as defined in the text and inserted into (gUHSJ, 
with p{m) taken from (|42|) in all three cases. 

for comparison. The three results for the r-values of the characteristic points A n and S n - 
cf. eq. (JHJ) for the relation between n and the winding number N - may be summarized as 
follows: 



where c n corresponds to a n , namely cq = l + mp(m), c\ = l+m/4, and c 2 = 1. In the regime 
of low m, all three choices of a give the same frequency ou(m). With p = po and the expansion 
of the log term we recover the result from the previous subsection: <jn = mr^ = lOjvv&Tjv 
with wn given by eq. (|28|) . For larger m the best agreement with the numerical integration 
of the original Lorenz equations is obtained from cq = 1 + mp(m): with the choice m = 0.75 
the prediction for the period of the pattern is 2P(mo) = 459.7. This is what was used to 
define the points S n and A n in Fig. namely, r N = 230 (N — \). 

If instead of N — j we take N — | in (jUJ), we obtain values r^{m) corresponding to the 
lower border of the r-bands. These are plotted in Fig.^J for iV from 1 to 9 and parameters m 
ranging from to approximately 1.25. The agreement with Fig. [2]is qualitatively convincing 
and quantitatively satisfactory. 




(44) 



with half period (in r) 




(45) 



24 



2000 
1500 
1000 



500 



500 1000 1500 2000 

FIG. 12: The curves (r N ,a N ) = ({N - \)P{m),mr N ) for < m < 1.25 and N = 1...9. The 
dashed lines are a = 0.75 r and a = r. 

V. BIFURCATIONS AND COEXISTENCE OF ATTRACTORS 

While the periodic orbits discussed in the previous two sections form the backbone of 
the structure in the (r, <r)-phase diagram on large scales of r and a, they do not represent 
the whole picture. In this section we report on the scheme of bifurcations that connects the 
different bands, and we also give a short discussion of the coexistence of attractors. 

A. Types of attractors: periods and symmetry 

In addition to the periodic attractors of periods 1 and 2, there are small bands of stability 
for periodic orbits of all periods, and of chaotic motion. This is summarized, using a color 
code, in Figs.^and|21 Within a given band, the attractors are the same type, i. e., they share 
the same period and symmetry property with respect to the reflection E, see (|5|). We use 
the notation "ps" for p-periodic symmetric attractors, and "pa" for p-periodic asymmetric 
attractors (whose mirror images under S are also attractors). 

As the structural features of the (r, er)-phase diagram are different in the regions of low 
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FIG. 13: (a, x)-bifurcation diagram at fixed r = 178, b = 8/3 for 4 < a < 65. The initial 
conditions are xq = yo = 0.001, Zq = r — 1 (red) and xo = yo = — 0.001, Zo = r — 1 (green). 
Attractors are asymmetric when they carry only one of the colors red or green. There are also 
symmetric attractors, given in black (superposition of red and green, as they are reached from 
both initial conditions). A few periodic orbits are shown in (y, z)-projection, of periods 2 in the 
upper right part, and periods 4 in the lower right; arrows indicate the corresponding branches of 
the bifurcation diagram. The conspicuous merging of the two major chaotic bands, as a is lowered, 
seems to occur at or very close to the value <r s (r) = Vbr = 21.8. 

and high a (separated roughly by cr s (r) = Vbr), we discuss them in turn. To begin with, we 
present in Fig. El a survey on the bifurcation behavior of attractors along the line r = 178 
for 4 < a < 65 . This cr-range bridges the two regions of small and large Prandtl numbers, 
respectively. It starts well below \fbr (= 21.8 for this r) and ends far above. The bifurcation 
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scheme along this line is shown in the left part of Fig. ED It records the x-values of the 
attractors in the Poincare surface of section z = r — 1. Red points were computed with initial 
conditions x = y = 0.001, z = r — 1, green points with x = y = — 0. 001, z = r — 1. 
When the colors do not mix, they represent asymmetric attractors; symmetric attractors are 
reached from both initial conditions and appear in color black. Starting with the highest a, 
we observe a small black segment of symmetric orbits with period 2, i.e., the attractors 
are type 2s. The (y, z)-projection for a = 63.5 is shown at the right side. Lowering a to 
approximately 62, the 2s-orbit splits into two 2a-orbits, one red, the other green, mirror 
images of each other under S. This is illustrated for a = 57. Each of these two attractors 
undergoes a period doubling scenario, 2a — > 4a — > 8a — > . . . , corresponding to red, green, 
blue bands in Fig. [TJ until chaos appears at about a = 50.5. When a is further decreased, 
the inverse cascade of band merging as described in ^(J is observed. It ends at a ~ 48.5, 
where the two bands of asymmetric chaos merge into one band of symmetric chaos, thereby 
restoring the symmetry under E. 

As a falls below approximately 41.2, a period 4 window emerges by an intermittency 
scenario. At first the attractor is symmetric, of type 4s, as shown for a = 40.7. Again, this 
orbit undergoes a pitchfork bifurcation into two asymmetric orbits of period 4, shown in red 
and green for a = 39.7, and then the period doubling sequence takes over for both of them. 
In Fig. E this window appears as the dominant green band in the white chaotic region. 

Fig. El shows the upper part of the bifurcation diagram in Fig. ED in greater detail, for 
the intervals 45 < a < 55 and 8.0 < o < 9.5. Except for their different scales the pictures 
are much like mirror images of each other. In both cases the scenarios coincide in many 
aspects (though not completely) with the bifurcation scheme of the antisymmetric cubic 
map 

x n+1 = f(x n ) = (1 - c)x n + cxl , -1<x<1, (46) 

shown in Fig. i n the ranges 3.2 < c < 3.4 and 0.25 < x < 0.8. This similarity was 
pointed out and explored using symbolic dynamics in [t| and jj. The reason for the good 
correspondence seems to be that (|46|) is the simplest polynomial 1-dim map that shares 
with Q a reflection symmetry. The invariance with respect to £ is here represented by the 
symmetry under x t— > Y,x = —x, namely S _1 /S = /. This modifies the standard period 
doubling scenario of the quadratic map x i— > ax(l — x) in such a way that the first bifurcation 
of a symmetric orbit of period p is not a period doubling ps — > (2p)s but rather a pitchfork 
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FIG. 14: Blowups of two pieces of the upper part of Fig. 1131 Left: 8 < a < 9.5, upper left branch; 
right: 45 < a < 55, upper right branch. The two pictures are like mirror images of each other. 
Again, red and green indicate the asymmetric attractors, approached from the two S-images of 
the initial values. 

bifurcation where ps breaks up into pa and its mirror image E(pa). Only then does period 
doubling start: pa — > (2p)a — > (4p)a 
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FIG. 15: Detail of the (c, x)-bifurcation diagram of the cubic iteration (j46fl . 

The character of the bifurcation diagram Fig. EH stays the same as we increase r and take 
any a-segment that connects the lowest two 2s-bands across the line <J s (r) = y/br. 
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Let us now consider what kind of regularity we can detect in the region of large r and a, 
in the onion like structure of Fig. |21 A representative picture of the bifurcation behavior is 
obtained by considering any parameter interval that connects two consecutive bands of equal 
symmetry. We take the line from Ai to A 2 and show the corresponding (r, ^-bifurcation 
diagram in the upper left part of Fig. El with the two initial conditions Xq = y$ = 10~ 5 
(red) and xq = yo = — 1CT 5 (green); zq = r + 1000 in both cases. (This large value of z 
is to a large extent arbitrary; however, it was observed that chaos coexists with periodic 
orbits, and that initial conditions with z of order r or lower tend to be captured by the 
latter.) The Poincare surface of section is zq = r — 1 as usual. Only the branch with 
x > is shown: the complete picture contains the mirror image x h- > —x, with colors red 
and blue interchanged. The small windows in Fig. El present a number of characteristic 
attractors in (y, z) -project ion, arrows connecting to the respective points on the bifurcation 
diagram. On the upper branch starting from A 2 , we have asymmetric periodic attractors 
of type la. These undergo a period doubling scenario (orbit 2a, etc.), including the inverse 
cascade (with asymmetric chaotic orbits), before they merge with their image under E and 
yield two symmetric chaotic orbits. This scenario is illustrated with the orbits in the upper 
right part of the picture. Upon lowering r, the chaotic attractor suddenly disappears, and 
an orbit of type 2s emerges, characteristic of the red band with center S 2 . This symmetric 
orbit splits into two asymmetric orbits of type 2a, in a pitchfork bifurcation like in Fig. fH?l 
(whose rightmost part belongs to the band Si rather than S 2 ). Next, period doubling sets in 
again. The chaos in the range of inverse cascading is at first asymmetric (though different 
from the one in the yellow band), but turns into symmetric chaos at the end point of the 
band merging. Finally, an attractor of type la emerges from this chaos and gives rise to a 
picture around Ai similar to the one with which we started at A 2 . 

We have not investigated the mechanism by which the "new" stable periodic orbits are 
generated when the chaotic motion becomes unstable upon lowering r. In fact, they coexist 
with the chaotic attractors before they are seen in Fig. El We cannot tell whether they 
emerge in an intermittency scenario. It may be that their range of existence, ignoring the 
question of their stability, extends all the way up to infinity. 

Fig. El is an extended version of the bifurcation diagram in Fig. El The range of r is 
from 50 to 2050, and instead of the ordinate x we plot y/r, again for z — r — 1. Remember 
that for a given periodic orbit, y reaches its maximum value when z is near r. 
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FIG. 16: Upper left: (r, x)-bifurcation diagram with m = a jr = 0.75 and r = 400 . . . 870, b = 8/3. 
The initial conditions are zq = r + 1000 and xq = yo = 10~ 5 (red) or xq = yo = — 10 -5 (green). The 
8 small windows show attractors corresponding to the parameter values identified by the arrows, in 
(y, z)-projection. Following the sequence in clockwise direction, starting at the top, the bifurcation 
scheme between A2 (r = 863) and Ai (r = 403) is illustrated: la-orbits (yellow band with center 
A2, see Fig. |2J) undergo period doubling, then become chaotic. Out of this chaos, 2s-orbits emerge 
(red band with center S2) and undergo a pitchfork bifurcation into 2a and its E-image; these 
bifurcate according to the period doubling scenario before chaos sets in, first asymmetric, then 
symmetric. Finally, la-orbits (of the Ai-band) emerge and the scenario starts anew. 

The picture contains two messages: (i) the repetitive nature of the pattern in Fig. |21 
extends to the bifurcation scenarios, and (ii) the value of y/r at its maximum, for periodic 
attractors, seems to tend to a constant p = p(m)/m as m = const and r — > 00. This was an 
essential assumption in the scaling scheme that we used in the analysis of Section TlVBI The 
value p = 0.37 that our theoretical investigation gave for m = 0.75 is close to the observed 
limit. 
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FIG. 17: (r, ?//r)-bifurcation diagram with m = a/r = 0.75 and 50 < r < 2050. The picture shows 
that the repetitive structure of Fig. [2] includes the bifurcation schemes. 



B. Coexistence of attractors 

It has already been mentioned several times that different attractors may coexist. For 
example, whenever an attractor A is not in itself symmetric under £, it coexists with its 
partner EA. This applies in particular to the fixed point attractors P± in the black region, or 
to the la-attractors in the yellow bands. Consequently, there must be at least four coexisting 
attractors where the yellow bands extend into the black region, as can be seen in Fig. |2j 
A less trivial case is analyzed in Fig. The parameters (r, a) = (108.82, 3.6) are located 
in the lower part of Fig. Q where the green asymmetric period 4 band overlaps with the 
black region. The left part of Fig. [TSl shows the Poincare surface of section z — r — 1, with 
initial conditions (xq, yo) colored pixel by pixel according to the final state to which they are 
attracted: white or blue if the final state is the fixed point P± = (±16.96, ±16.96, 107.82), 
respectively, red or gray if it is one of the two 4a-orbits. No other attractor besides these 
four was found. The blue and white basins of the point attractors P± appear to have smooth 
boundaries whereas the boundary between the basins of attraction of the two 4a-orbits is 
highly entangled. 



31 



FIG. 18: Left: Basins of attraction of 4 coexisting attractors, for (r, cr) = (108.82,3.6). Initial 
conditions in the (xo, yo)-plane, —91 < xq,uq < 91, with zq = r — 1, are scanned pixel by pixel 
and given a color depending on the attractor to which they are asymptotic: white for P + , blue for 
P_, red or gray for the asymmetric period 4 orbits shown on the right (abscissa y and ordinate z). 
The fixed point attractors P± are marked as black dots. All intersection points (8 each) of the two 
period 4 attractors with the Poincare surface z = r — 1 are emphasized. 

Another case of coexisting attractors is documented in Fig. EH where the parameters 
(r, a) = (1200,255.64) are from the chaotic band at the upper rim of the lowest yellow 
band in Fig. |21 Fig. demonstrates how a symmetric chaotic attractor coexists with two 
attractive asymmetric period 1 orbits. No other attractors have been detected. The basin 
boundaries seem to be complicated but not fractal; they are spiralling out from the unstable 
fixed points P±. Each point is colored according to which attractor absorbs it. Each basin 
trivially contains its respective attractor. The location of their intersections with the plane 
z = r — 1 is marked by black dots. 
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FIG. 19: Left: Basins of attraction of three coexisting attractors for (r,a) = (1200,255.64). Initial 
conditions with —355 < xq,uq < 355, zq = r — 1 are scanned and colored according to their final 
fate. The three attractors are shown on the right, in (y, z)-projection: two asymmetric orbits of 
period 1 (with the blue and green basins of attraction), and one symmetric chaotic attractor (in 
black, with the white basin of attraction). The intersections of the attractors with the Poincare 
surface are emphasized by arrows. 

Inspection of Figs. and EH suggests that a scan of initial conditions along the line 
Xq = yo, z = r — 1, might catch all (co-)existing attractors. Thus instead of scanning a 
two-dimensional set of initial conditions, we scan along this diagonal line only and use the 
second dimension to superimpose these scans for a range of parameters a. Such a "cx-stack" 
is shown in Fig. 1201 for r = 1200. The abscissa is xq = yo in the range from to 120, and 
the ordinate is a from to 400. Comparison with Fig. |2 shows that this a-range covers the 
bands containing the points Si, A 1; and S2, in addition to the low-u structure with regular 
behavior at the very bottom and chaos in a range around \/br = 56.6. 

The coloring scheme is somewhat arbitrary and was designed so that the coloring of 
different attractors could be done automatically. For a fixed value of a different color 
designates different attractors, i.e., the number of different colors on a scan-line is the 
number of coexisting attractors. Same color for adjacent points with different a implies 
that the corresponding attractors are deformations of each other. Color white is reserved 
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for chaotic motion. 

Let us concentrate on the a-ranges which are marked as Sf a , Sf s , Ci, etc. Sf s is the 
a-interval (126.5, 154.8). The attractors are symmetric orbits of period 2 and type Si, see 
Fig. and there are no other attractors; hence there is only one color. For smaller a in the 
range Sf a , 113 < o < 126.5, the attractors have undergone a pitchfork bifurcation and exist 
in pairs of asymmetric orbits of period 2. Thus we see two colors, corresponding to the two 
basins of attraction. At larger a, on the other hand, in the range Ci = (154.9, 190.3), the 
attractors of type 2s coexist with chaotic attractors; again we observe two colors. 

The band with orbits of type Ai exists in the range A\ a = (207.3, 250.8). As the attractors 
are asymmetric orbits, they occur in pairs, and two colors appear (black and magenta). At 
lower a, we have a range A\ a = (190.3, 207.3) where period doubling has increased the period 
but not the number of attractors; thus there are again two colors. At larger a, in the range 
C2 = (250.8,283), the two asymmetric periodic attractors coexist with a chaotic attractor, 
and three colors are needed. Fig. is an example from this range. 

If a is further increased, we enter the regime which is discussed in Fig. E3 The bifurcation 
tree of the S2-band leaves a complicated trace in Fig. 1213 but then there comes a range where 
a 2s-orbit is the system's only attractor. From thereon, the pattern seems to repeat itself, 
but the details on the x = y axis become more and more intricate. 
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FIG. 20: o"-stack of domains of attraction, for r = 1200. For a in the range between and 400 
(ordinate) initial points along the diagonal xq = yo (and in the Poincare plane z = r — 1) are 
given a color according to which orbit attracts them. White, in particular, means there are chaotic 
orbits. The values xq = yo are taken from the range (0,400). The features of the cr-ranges Sf a , 
Sl s , Ci, etc. are explained in the text. 
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